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1 .  BACKGROUND 


In  order  to  characterize  an  unknown  object  by  remote  sensing,  certain 
electrical  properties  of  the  object  may  be  used  [1-9] •  In  particular, 
the  impulse  response,  ramp  response  and  the  natural  resonant  frequencies 
of  the  object  have  been  proposed.  In  order  to  accurately  probe  the 
object,  one  must  illuminate  it  with  a  band  of  frequencies  whose  wave¬ 
lengths  are  approximately  the  same  dimensions  as  the  overall  length 
of  the  object.  Since,  in  general,  the  size  of  the  object  is  not  known 
beforehand,  it  is  necessary  to  illuminate  the  object  by  a  broadband 
signal  such  as  a  waveform  which  is  an  approximation  to  an  impulse. 

The  transmitted  pulse  induces  electrical  currents  on  the  object  and 
reradiates  the  Incidental  energy.  The  problem  Is  then  to  obtain  the 
electrical  properties  of  the  object  from  the  transmitted  and  received 
time  domain  waveforms.  This  area  has  been  investigated  by  the  researchers 
at  the  Ohio  State  University  [1-5] . 

In  recent  years  with  the  advent  of  the  singularity  expansion  method, 
efforts  have  been  directed  towards  characterizing  an  object  by  its 
natural  frequencies  -  poles  and  zeros  [6-8],  For  this  purpose,  Prony's 
method  has  been  applied  with  much  success  in  analyzing  measured  impulse 

responses  with  high  signal  to  noise  ratio.  However,  it  is  difficult  to 
extend  Prony's  method  to  arbitrary  Input  and  arbitrary  output  wave¬ 
forms.  Also  the  SEM  representation  does  not  account  for  the  impulsive 
portton  of  the  system  time  behavior.  Hence,  the  complete  Impulse  response 
contains  more  Information  than  the  SEM  poles.  One  of  the  objectives  of 
this  paper  Is  to  obtain  the  Impulse  response  In  the  time  domain  when 
arbitrary  input  and  output  waveforms  are  given. 

More  recently,  the  pencil  of  function  method  [10-14]  has  been  applied 
with  much  success  to  obtain  the  natural  frequencies  of  objects  from 
measured  arbitrary  Input  and  output  time  waveforms.  Both  Prony's  and 
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the  pencil  of  function  method  are  frequency  domain  techniques.  The 
accuracies  of  determining  Impulse  responses  by  either  of  these  techniques 
depend  directly  on  the  time  domain  data  available.  Also  a  significant 
problem  with  any  such  parametric  scheme  is  the  need  to  estimate  the 
number  of  dominant  poles  and  zeros.  A  poor  estimate  can  lead  to  large 
errors  in  the  results. 

Finally  as  discussed  In  the  next  section,  it  is  not  easy  to  obtain 
the  impulse  response  from  given  input  and  output  time  domain  waveforms 
using  Fast  Fourier  Transform  techniques  t9l. 

The  objective  of  this  paper  is  to  present  two  techniques  for  obtaining 
the  impulse  response  of  an  object  directly  in  the  time  domain  from 
measured  input  and  output  waveforms. 

2.  INTRODUCTION 

Let  x(t)  be  the  input  to  a  time  invariant  causal  linear  system  which 
is  characterized  by  its  impulse  response  h(t).  The  corresponding  output 
y(t)  is  given  by 

CO 

y(t)  -  /  x(t-T)h(t)  dr  (I) 

o 

One  approach  for  solving  (1)  for  the  impulse  response  involves  the 
Laplace  Transform.  The  result  Is 

Y(s)  -  H(s)«X(s)  (2) 

where  X(s),  Y(s)  and  H(s)  are  the  Laplace  Transforms  of  x(t) ,  y(t)  and 
h(t),  respectively.  The  impulse  response  h(t)  is  obtained  by  taking  the 
Inverse  Laplace  Transform  of  (2)  to  yield 

h(t)  *  2*J  (,  ?{sf  exp(st)  ds  0> 
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where  T  Is  the  Bromwich  contour.  Numerical  problems  are  frequently 
encountered  with  this  approach  when  dealing  with  real  data  [9].  For 
example,  when  using  Fast  Fourier  Transform  techniques,  the  division 
blows  up  If  X(jw)»0  In  the  pass  band. 

Two  time  domain  methods  are  presented  in  this  paper.  One  Is  the 
method  of  synthetic  division  while  the  other  is  the  method  of  least 
squares.  The  former  is  computationally  straightforward  whereas  the 
latter  Is  computationally  more  stable. 

in  this  presentation  y(t)  and  x(t)  are  assumed  to  be  discrete  time 
signals  with  identical  sampling  rates.  The  problem  is  to  process  the 
discrete  data  points  so  as  to  estimate  the  corresponding  discrete 
impulse  response  without  introducing  additional  information  like  the 
augmented  high  frequency  response  [9]. 

Let  the  tth  samples  of  y(t),  x(t)  and  h(t)  be  denoted  by  y(,  x{, 
and  hJf  respectively.  Also,  assume  that  y(t)  and  x(t)  are  available 
as  the  finite  time  series  {y)  -  {yQ,  y^,  Y2’,,,,yN^  and  M  *  {x0»  xi * 
x2, • • • The  object  Is  to  estimate  the  finite  Impulse  response 
time  series  {h}  •  {h^,  h^,  h2,...,h^}. 

For  convenience,  the  following  assumptions  are  made: 

(a)  h,  -  0  for  I  >  L. 

(b)  Xj  Is  observed  for  0  J5.  J  M  and  Is  not  Identically  zero. 

(c)  Xj  -  0  for  0  >  J  and  J  >  H. 

(d)  yk  Is  observed  for  0  jc  k  £  N  (-M+L) 


3.  BA  HU'S  TECHNIQUE  (SYNTHETIC  DIVISION)  [15-16] 

This  method  Is  applicable  for  all  positive  Integer  choices  of  N, 

M  and  L.  The  sequence  corresponding  to  h(t)  Is  written  symbolically  as 


fh}  _  _  {vo,yl’y2’“-*vN}  r  .  n 

{h}  itrm  rx0,x,,x2 . xmj  “°vhi . hl 


(4) 
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where  the  division  of  sequence  {y}  by  sequence  {x)  is  carried  out  by 
synthetic  division. 

Justification  of  the  synthetic  division  operation  is  obtained  from 

the  theory  of  2-transforms.  The  input  and  output  z-transforms  can  be 
written  as 


-1  -7  -N 

♦y  Z  +...+yyz 

4  N 

(5) 

:  .  .+x^z 

(6) 

01 

h(z)  is  then  the  ratio  of  the  two  polynomials  given  by  (5)  and  (6).  Since 
two  z-polynomiais  are  divided  in  the  same  manner  as  algebraic  polynomials, 
the  operation  of  synthetic  division  is  justified. 

Alternatively,  a  discrete  approximation  to  (I)  is 

y At  |  h,  ,x,  (7) 

J  i-0  J 


i 


where  At  is  the  sampling  interval.  A  computationally  compact  form  for 


hj  is  then  given  by 


h,  * - a* 

j  xnAT 


j-ATj,xiYi 


(8) 


where  it  is  assumed  Xgi*0.  The  values  for  hj  obtained  in  this  way  differ 
from  those  obtained  by  synthetic  division  by  the  constant  factor  1/At. 

Note  that  both  of  these  methods  are  applicable  to  discontinuous 
inputs  and  outputs. 

Example  1;  As  a  simple  example,  consider  the  rectangular  pulse  {x}  ■ 

{1,  1,  1}  to  be  an  input  to  a  system  whose  output  is  the  triangular  pulse 
{y }  -  {1 ,2,3,2, l}.  Using  synthetic  division,  as  implied  by  (4)  the 
unknown  transfer  function  {h}  is  given  by 


hU)  -  i-2 

1+2  '+2 


(9) 


This  implies  {h}  »  0,1,1). 

Alternatively,  if  it  is  assumed  At  •  1,  application  of  (8)  results 


In  {h>  •  0,1,1).  This  agrees  with  the  solution  obtained  by  synthetic 
division. 


k.  EFFECT  OF  MEASUREMENT  ERRORS 


The  effect  of  measurement  errors  on  the  synthetic  division  scheme  is 
now  considered.  Assume  that  the  input  measurement  is  corrupted  by 
additive  noise  p(t)  so  that  the  total  recorded  input  is  x+p.  Also 
assume  that  the  output  measurement  is  corrupted  by  noise  q(t)  such  that 
the  total  measured  output  is  y+q.  Using  samples  of  the  corrupted  input 
and  output  measurements,  an  estimate  of  the  impulse  response,  denoted 
by  Ti(z) ,  is  given  by 


fi(z) 


y(z)+q(z) 

TO157 


(y0+v  ,+y2z  2+-->+Ynz  N)+(q0+qi2  1+---+(iN2  N) 

.1  .  |  -  — |y|  _ 

(xq+XjZ  +...+xmz  )+(p0+PjZ  +...+pmz  ) 


(10) 


in  this  section  the  relationship  of  the  estimate  {h}  to  the  true 
impulse  response  {h}  is  explored  where  it  is  assumed  that  the  measure¬ 
ment  errors  are  small  and  that  the  elements  of  {q>  and  {p}  which  pertain 
to  two  different  measurements,  are  statistically  independent  zero  mean 
random  variables. 

The  expression  in  (10)  can  be  rewritten  as 


h(z)  -44-^41 


x  (z) 


For  those  values  of  z  for  which 


(11) 


<  1,  H  + 


may  be  expanded  in  a  Binomial  series  to  yield 

[1  + 


w;  w;  i- . , .  . ,  j&f )» 


(12) 


Xq+X^Z  ♦. .,+x^z 


It 


follows  for  these  values  of  z  that 


(13) 
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The  expected  value  of  a  random  sequence  is  defined  to  be  the  sequence 
whose  elements  are  expectations  of  the  original  elements.  Since  each 
element  represents  a  sample  value  of  the  waveform,  the  expectation  of 
the  element  corresponds  to  an  ensemble  average  of  the  waveform  at  the 
Instant  of  the  sample  value.  Then  the  expected  value  of  the  estimate 
h(z)  is  given  by 


£[q(z)]  -  £’[q()+q1z”1+q2z*2+...3 

-  Flqgl+Elq^z  1+£,[q2]2"2+... 

-  0  (15) 


where  use  has  been  made  of  the  assumptions  of  zero  mean  and  statistical 
Independence  of  the  random  variables  P|  and  q^.  Similarly, 

E  Ep(z) ]  -  0  (16) 

and 

E[p(z)l2  -  EIpq^PqPjZ  ^p^z  2+... 

2  -2  -3  -L  _r 

+p,  z  +2P,P2*  ^+...+P2z  +2p2P3z 

2  2  -2  2  -ii 

*  o  *  +  a  z  +  a  z  H+  ...  (17) 

P0  Pi  P2 

where 

°p,2  ’  *Ipil  *  (,8) 

Consequently  0  2  +  0  2^2  + 

Elhi z)]  -  h(z)«[l+  - -  ] 

(Xq+XjZ  +. ..+x^z  ) 

-  h'0  +  hjz'1  ♦  h£z'2  +  ...  (19) 
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where  h(z)  ■  y(z)/x(z),  Thus  If  the  variances  of  the  Input  measurement 
errors  are  known,  approximate  values  for  the  expectation  of  the  estimate 
of  the  time  domain  Impulse  response  fi(z)  can  be  evaluated.  Notice  that 


If  pj  =0  for  all  I  (l.e.  the  Input  measurement  Is  free  from  any 
additive  noise)  then  the  estimate  Is  unbiased  even  though  a  measurement 
error  Is  still  associated  with  the  output.  However  If  Pj  t  0  then  there 
Is  a  bias  In  the  solution  which  Is  given  as 


E[h(z)]  -  h(z)  -  h (z) 


r  \</**  --  \ 

\  (x0+x1z“,+...+xmz'M)2  J 


By  definition,  the  variance  of  a  sequence  is  the  sequence  whose 
elements  are  the  variance  of  the  original  elements.  Each  variance  is 
the  variance  of  the  original  random  waveform  at  the  corresponding 
sampling  instant.  Note  that 

Mz,  -  IPIU)]  -  »(,)  I  *{f&}  2-  ^  1 


where  a  (z) 
P 


♦  ...  higher  order  terms 

2  2  -2  2  -k 

a  +az  +oz  +... 


(20) 


'0  n  r2 

As  a  first  approximation 

K(z)  -  Flfi(z)]  *  q(-z)  V(V)(?if— 


-1  -2 

co  +  c\2  *  v  + 


(21) 


Note  that  (21)  corresponds  to  the  zero  mean  sequence  {c^.Cj  .c^,...}. 

It  follows  that  the  variance  of  this  sequence,  which  also  equals  the 

variance  of  h(z),  results  in 

o2(fi(z)]  -  Slc^l+fflc^Jz  l+E[c22]z  2+... 

2  2-1  2-2 
■  ffj  ♦  Oj  Z  ♦  Oj  Z  +  ... 

2  2 

where  Oj  »  ff[cj  ). 


(22) 


Use  of  (22)  and  (16),  in  conjunction  with  Tshebyshev's  inequality 


t 1 7] »  yields  the  inequality 

P[hj  -  kOj  <  fij<  hj  *  kOj]  _>  1 - (23) 

k 

where  P [ . )  denotes  the  probability  of  the  argument  in  the  brackets  and 

Kj  is  the  element  In  {fi}.  Equation  (23)  is  valid  for  small  measurement  errors . 

In  summary,  if  the  first  and  second  order  statistic  of  the  random 

sequences  {p}  and  (q)  are  known,  then  (22)  and  (16)  yield  approximate 

expressions  for  the  variance  and  mean  of  the  sample  values  of  the 

extimate  of  the  time  domain  impulse  response  obtained  from  (10). 

Example  2:  Consider  a  system  with  impulse  response  (h)  ■  {1,1}  for 

which  the  input  is  {x}  *  {1,1}  and  the  output  is  {y}  *  {1,2,1}.  Assume, 

for  the  purpose  of  illustration,  that  the  input  and  output  are  each 

contaminated  by  zero  mean  additive  measurement  errors  which  take  on 

the  values  +0.01  or  -0.01  with  equal  probability.  Note  that  the 

2  -A 

measurement  errors  have  zero  means  and  variances  o  *  10  .  Consider 

a  particular  measurement  error  which  results  in  the  observations 
{x}  -  {1.01,0.99}  and  {y}  -  {0.99,  2.01,  1.0}.  In  this  example,  the 
effects  of  the  measurement  errors  on  the  extimate  {h}  are  explored. 

From  (10)  the  estimate  of  the  impulse  response  is  obtained  as 
~  -1-2 

R(z)  -  -  -99+2-°Vz  +* — »  . 9802+1. 0293z'1-.0190z“2+.0190z"3  (2 A) 

xU'  1 .01+.99z” 

{R}  -  {.9802,  1.0293,  -.0190,  .0190,  ...}  (25) 

The  measurement  errors  are  seen  to  cause  the  elements  in  {R}  to  differ 
from  those  in  {h}  and  to  cause  {R}  to  be  an  Infinite  sequence  whereas 
{h}  was  finite. 
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By  using  (24)  the  expected  value  of  the  estimate  is  found  to  be 

£[fi(z)J  -  (l+o2)+(l-a2)z“l+2aV2-2oV3+... 

-  1.0001  +  0.9999z_1  +  .0002z“2  -  .0002z"3+  ...  (26) 

It  is  seen  that  the  mean  value  of  the  estimate  is  biased.  This  is  due 
to  measurement  errors  in  the  input. 

The  variances  of  the  estimate  are  obtained  from  (21)  and  (22). 
Observe  that 

s(«wi  -  aULijMikto 

-  (q0-p0)*(q1-p|-q0)z‘l+(q2-q1+q0)2"2-(q2*q1-Ki0)z'3.  ...  (27) 

It  follows  that 

o2[fi(z)]  -  I^0"P0^2J+  £'I(cl|-Pj"cl0)2]z"1+Et(q2-q,+q0)2z‘2 

+  £[(q2-qj+q0)2]z"3+  ... 

-  o2[2+3z_1+3z”2+3z’3+...]  (28) 

2 

Consequently  the  variance  of  the  first  element  in  {K}  is  2 o  while  it 
2 

is  3 o  for  the  remaining  elements. 

Letting  k*3  in  (23),  Tshebyshev's  inequality  yields 
P[|fi0- 1.0001 |  <  3(.Ol4l)]  >  8/9 
P[|R, -0.9999|  <  3(. 0173)3  >  8/9 

P[|fi2-  .0002 1  <  3 ( - 0 1 73) 3  >  8/9  (29) 

Observe  that  the  estimates  in  (25)  are  within  the  3<?  of  the  expected 
values. 

A  disadvantage  with  the  synthetic  division  scheme  is  that  the  error 
tends  to  build  up  with  the  length  of  the  sequence. 
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5.  A  LEAST  SQUARES  APPROACH  TO  THE  ESTIMATION  PROBLEM 


If  X(s)  [which  is  the  Laplace  Transform  of  x(t)]  contains  a  zero 
In  the  right  half  plane  (RHP)  [  i.e.  non  minimum  phase],  a  serious 
computational  difficulty  may  arise  with  the  synthetic  division  scheme 
because  H(s)  may  then  contain  a  pole  in  the  RHP  [ 1 8] .  This  results  in 
numerical  instability.  With  a  stable  system  the  RHP  zero  of  X(s)  is 
cancelled  by  the  identical  RHP  in  Y(s).  However,  in  practice,  round-off 
and/or  truncation  errors  in  computing  X(s)  and  Y(s)  prevent  this 
cancellation.  The  synthetic  division  approach  does  not  then  converge. 
Should  this  occur,  an  alternate  approach  must  be  used.  This  section 
develops  the  least  squares  estimation  of  the  time  domain  impulse 
response  from  measured  input  and  output  data. 

With  reference  to  (7),  the  response  of  a  discrete  time  system  may 
be  written  as 


y 


J 


(30) 


The  convolution  in  (30)  may  be  expressed  in  matrix  form  as 


Y0 

xo 

0  0 

0 

...  0 

V 

Yi 

• 

A 

XJ 

X 

o 

o 

0 

...  0 

hi 

• 

• 

X2 

X  j  •  * 

• 

•  •  • 

• 

Yn 

NXl 

XM 

Vi  ’ 

• 

•  •  • 

.V 

-  rn 

0 

0  .  . 

• 

•  •  xM 

— 

n 

NXL 

[V’]  -  [XJ[H]  (31) 


where  Y'  and  H  are  NX1  and  LX1  column  matrices,  respectively,  and  X  is  an 
HXL  rectangular  matrix.  Equation  (31)  can  be  solved  for  [H]  by  first 
premultiplying  both  sides  by  [X]^  (l.e.  transpose  of  [X]).  Then 
[X]T[X] [H]  -  [X]T[Y'3  (32) 
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It  follows  that 

[H]  «  {[X]T[X]}‘,[X]T[Y']  (33) 

It  is  interesting  to  note  that  [X]^[x]  is  the  autocorrelation  matrix  of 

the  available  input  time  samples.  Thus,  [X]T[X]  is  a  symmetric  positive 
seml-def Inlte  matrix.  In  fact  [X]^[X]  is  a  Topittz  matrix  and 

therefore  Is  Invertible  by  means  of  a  recursive  procedure  requiring 

2  T 

L  operations.  The  Inversion  of  [X]  [X]  Is  numerically  highly  unstable 

as  has  been  pointed  out  by  Ekstrom  [191*  He  suggested  a  Singular  Value 

Decomposition  technique  for  the  inversion.  In  this  procedure  [X]^[X] 

is  converted  to  a  diagonal  matrix  by  a  similarity  transformation.  The 

eigenvalues  are  then  filtered  (i.e.  eigenvalues  below  a  certain  small 

J 

amount  are  set  to  zero).  The  whole  process  of  premultiplying  by  [X] 
and  then  obtaining  the  singular  value  decompos I t I on  is  highly  time 
consuming.  In  this  paper  an  alternate  approach  is  suggested  based  upon 

the  conjugate  gradient  method  [20]. 

With  the  conjugate  gradient  method  it  is  not  necessary  to  form  the 

matrix  [X]^[X].  The  matrix  equation  [X] [H]  ■  [Y1]  is  solved  by  the 
following  procedure: 

By  usingany  initial  guess  [H]q  for  the  solution,  the  matrices  [R]g 
and  [P]q  are  generated  where 


[R]0  -  [X][H]0-[Y'] 

(W 

and 

[p]0  -  “ [x]t[r] . 

The  (n+2)**1  estimate  is  then  given  by 

(35) 

lH]n-M  "  [HK  +  lnlpln 

(36) 

where 

|[XJTIRJ  |2 

t  -  - 5- 

"  IWIP]„|2 

(37) 
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(38) 


and  [R]  and  [P]  are  obtained  according  to  the  relations 
n  n 

W.*Wrl’V|WMr, 

and 

[Pln  -  " [x]T[R]  +  q„  .[p]n.  (39) 

n  n  n- 1  n  i 

where 

|[x]tIr]  |2 

q  «  “  - e—2 - 7  .  m 

IM 

The  advantage  of  the  conjugate  gradient  method  is  that  ft  Is  an 
iterative  scheme  which  usually  yields  excellent  results  within  J  Iterations 
where  J  Is  the  number  of  Independent  eigenvalues  of  [X]^[X].  It  has  been 
shown  that  the  conjugate  gradient  method  converges  to  a  solution  even 
when  [X]T[X]  Is  singular.  For  this  case  the  iteration  process  Is 
terminated  once  the  solution  stops  converging  (l.e.  begins  to  oscillate). 
This  method  has  been  used  successfully  In  the  area  of  Image  processing 
for  a  similar  deconvolution  problem  [21]  .  However,  an  analytical 
expression  is  available  for  the  total  number  of  Iterations  beyond  which 
the  Iterative  procedure  can  be  terminated  [25]. 

The  error  analysis  performed  in  the  previous  section  is  now 
repeated  in  order  to  find  mean  and  variance  of  the  solution  obtained  by 
the  least  squares  method  when  measurement  errors  are  introduced. 
Corresponding  to 
[X|  [H]  -  [Y'J 

mi 

[H] ,  an  estimate  of  [H] ,  satisfies  the  matrix  equation 

mi 

[X+P]  [H]  -  [Y'+Q]  (1,1) 

where  [P]  and  [Q]  account  for  the  measurement  errors  and  are  defined  as 


By  assuming  the  measurement  errors  are  small,  an  estimate  of  the 
impulse  response  is  obtained  as 

[H]  -  {[X+P]T[X+P]f '[X+PlV'+Q] 

-  {[X]T[X]}",{[I]  -  ([X]T[Xjf ’{A]}  [X+PjV+Q] 

+  higher  order  terms  (A4) 


where 


II]  «  identity  matrix,  and 

[A]  -  [X]T[P]  +  [P]T[X]  +  [P]T[P]. 

The  expected  value  of  the  time  domain  impulse  response  is  given  by 

£[H]  =  {[X]T[X]}’,{[I]-ax]TtX])"1[aJ]}[X]T[Y'] 

where  [o2]  -  E{[P]^[P]}  .  (*»5) 

P 

Also  note  that  when  measurement  errors  are  associated  with  the 
input,  then  there  is  a  bias  in  the  solution  which  Is  given  by 
H  -  Elifi  S  {[XlT[X]r2t0p][X]TtY'}. 

Similarly  the  variance  of  [H]  is  defined  to  be  a  column  vector  whose 


elements  are  the  variances  of  each  of  the  elements  in  [H].  Note  that 


[fl]  -  FIB]  -  {[xftX]}”1  [P]T[V]  +  [X]T[Q]  -  ([X]T[X])'1 

*{[X]T[P]  +  [P]T[X]>  [X]T[V] 

+  higher  order  terms 

•  [0]  +  higher  order  terms  (46) 


Then 

o2[H]  =  £ 


(47) 


From  (44)  an  estimate  of  the  time  domain  impulse  response  can  be  obtained 
at  each  instance  through  the  column  matrix  [R] .  The  expected  value  and 
the  variance  are  obtained  from  (45)  and  (47),  respectively,  through  the 
elements  of  the  column  matrices  E[R]  and  o  [R] .  Using  Tshebyshev's 
inequality,  as  in  the  previous  section,  the  probability  of  a  particular 
element  [R  J  of  the  matrix  [H]  lying  between  E[RjJ  ±  ka [fl ^ ]  is  greater 

than  (1 - jr-  ) ,  i  .e. 


PtllHj]  -  FlBjH  <  kotfj]}  >  i  -  -i-  . 

Example  3*  Consider,  once  again,  the  same  problem  stated 
From  the  problem  statement,  (41)  becomes 


"l.Ol 

0 

"o 

.99 

0.99 

0 

1.01 

0.99 

« 

X 

■t 

2.01 

1.0 

After  premultiplying  both  sides  of  the  equations  by  [X+P]T 
[X+P)T[X+P],  the  estimate  [H]  becomes 


(48) 

in  example  2. 


(49) 

and  inverting 


(50) 
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This  compares  to 


(51) 


by  the  synthetic  division  method. 

in  addition,  the  least  squares  approach  yields 


e[ h]  =  {[xftxir1  ([!)  -  {[x]T[xl}-,[cjp]) [X]T[Y'] 


-3pq  +2qQ  +qQ  -  q2 
“3Pj  ‘  +2q2 


Thus  the  variance  of  the  estimate  is 


(52) 


o2[B] 


1.6667 

1.6667 


(53) 


It  is  seen  that  (50)  is  within  the  bounds  dictated  by  (48)  for  the 
special  case  in  which  k*3;  l.e.,  the  estimate  is  within  the  range  given 


By  comparing  (50)  to  (25)  we  find  that  the  least  squares  method  yields 
a  better  estimate  than  the  synthetic  division  technique  for  this 
problem.  Also,  the  variance  of  the  least  squares  method  is  much  smaller 
than  the  variance  of  the  synthetic  division  method.  However,  the 
synthetic  division  scheme  is  computationally  much  simpler  than  the  least 
squares  method  and  yields  excellent  results  when  a  of  the  measurement 
error  is  very  small. 


15 


7.  OPTIMUM  INPUTS  FOR  WHITE  NOISE 


In  most  practical  cases,  the  Input  Is  known  with  great  accuracy. 
The  significant  errors  are  then  associated  with  measurement  of  the 
output.  The  question  then  arises,  "Given  p(z)=0  and  q(z)i<0,  which 
Input  x(z)  minimizes  the  variances  of  h(z)  In  (22)  and  ( 47)?' ' 

It  Is  apparent  from  (21),  (22)  and  (47)  that  the  variances  of 
h(z)  are  reduced  as  the  amplitude  of  x(z)  is  Increased.  In  a 
practical  situation,  however,  the  amplitude  of  the  Input  may  be 
limited  In  order  not  to  excite  system  nonl Ineartties.  Also,  the  Input 
amplitude  may  be  constrained  due  to  power  limitations.  Therefore  the 
above  question  Is  pursued  subject  to  a  constraint  on  the  Input  ampli¬ 
tude.  For  convenience,  the  constraint  is  chosen  to  be  on  the  mean 
square  value  of  x(z),  i.e. 

<|>xx(0)  -  autocorrelation  of  the  Input  x  of  lag  zero 


M 

'  M+T  2  xm 
rm«0 

Assuming  p(z)  *  0  and  q(z)  to  be  a  white  noise  sequence,  (47) 
becomes 

*%  T  _1 

(56) 


ct2[H]  -  crq2{[X]T(X]>"1 . 

Observe  that  [X]T[X]  Is  a  symmetric  positive  definite  matrix  with  each 


element  along  the  principal  diagonal  having  the  value  (M+l)$xx(0). 
Using  Levin's  results  [24],  the  variance  of  (56)  is  minimized  If  and 
only  If 


[X]TIX]  -  (M+U^tO)  •  [I] 

where  [I]  Is  the  Identity  matrix.  Hence,  provided  the  Input  Is 
deterministic,  It  follows  that  the  solution  for  x(z)  satisfies 


(58) 


♦„<0>  "  ° 

$  (£)  ■  0  for  0  <  £  <  L 

xx  — 

Hence,  the  Input  sequence  x(z)  which  minimizes  the  variance  is  white 
over  a  range  of  L  samples. 

8.  MINIMUM  VARIANCE  ESTIHATE  OF  THE  IMPULSE  RESPONSE  FOR  COLORED  NOISE 
For  the  case  in  which  the  noise  q(z)  Is  not  white,  as  was  assumed 
In  the  preceding  section,  then  Laplace  [22],  Gauss  [23]  and  Markov  [24] 
have  shown  that  the  minimum  variance  unbiased  estimate  for  the  Impulse 
response  Is  given  by 

IX]T[$  l*1  [X] [Hm1  -  tXlVl'V]  (59) 

qq  M  qq 

Observe  the  similarity  to  the  solution  given  by  (32).  In  (59)  the 
matrix  14>^]  represents  the  covariance  matrix  of  the  noise  sequence  q(z). 
The  subscript  M  In  [H  ]  Implies  they  are  Markov  estimates. 

n 

Equation  (6l)  can  be  solved  In  a  very  computationally  efficient  way. 
Since  (♦qq)  Is  a  non-singular  symmetric  covariance  matrix,  Its  Inverse 
can  always  be  factored  Into  the  following  form 

[♦«,]“ 1  -  [V]T[W]  (60) 

By  application  of  (60),  (59)  can  then  be  written  as 
[C]T[C][HmJ  -  [C]T[W] [Y* )  (61) 

where 

ICJ  £  [V][XJ  (62) 

Equation  (62)  can  now  be  solved  by  the  conjugate  gradient  method  as 
outlined  tn  section  5. 

If  the  noise  Is  white  then  (59)  Is  equal  to  (32)  and  the  solutions 
are  Identical.  If  the  noise  q(z)  Is  not  white,  the  Markov  estimates 
[H^]  are  much  more  difficult  to  compute.  However  the  disadvantages  of 
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dotng  more  computation  are  offset  by  the  fact  that  the  Harkov  estimates 

are  the  minimum  variance  unbiased  linear  estimates  of  [H]  (t.e,  linear 

In  [Y']).  If  the  noise  sequence  q(z)  Is  Gaussian,  It  can  be  shown  that 

[HJ  Is  an  unbiased  estimate  and  the  Cramer-Rao  inequality  holds  with 

n 

the  equality  sign  In  [25]. 

9.  CONCLUSION 

Two  techniques  have  been  presented  to  compute  the  time  domain  Impulse 
response  from  measured  time  domain  Input  and  output  waveforms.  The 
computation  Is  carried  out  directly  in  the  time  domain  without  the 
need  to  perform  any  Laplace  Transforms.  Error  bounds  are  presented  when 
the  measured  waveforms  are  contaminated  with  noise  and  the  statistics 
of  the  measurement  errors  are  known.  The  Input  required  under  an 
amplitude  constraint  to  minimize  the  variance  of  the  Impulse  response 
estimate  Is  also  discussed.  Finally,  a  generalized  least  squares 
technique  Is  presented  which  yields  a  minimum  variance  unbiased 
estimate  for  the  Impulse  response. 

Vork  Is  currently  under  way  dealing  with  numerical  examples  from 
measured  waveforms. 
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